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We present the first comparison between numerical relativity (NR) simulations of an eccentric 
binary black hole system with corresponding post-Newtonian (PN) results. We evolve an equal- 
mass, non-spinning configuration with an initial eccentricity e w 0.1 for 21 gravitational wave cycles 
before merger, and find agreement in the gravitational wave phase with an adiabatic eccentric 
PN model with 2 PN radiation reaction within 0.1 radians for 10 cycles. The NR and PN phase 
difference grows to 0.7 radians by 5 cycles before merger. We find that these results can be obtained 
by expanding the eccentric PN expressions in terms of the frequency-related variable x = (cj M) 2 ^ 3 
with M the total mass of the binary. When using instead the mean motion n = 2n/P, where P is 
the orbital period, the comparison leads to significant disagreements with NR. 



I. INTRODUCTION 

Tremendous progress towards detecting gravitational 
waves is being made by observational efforts such as 
LIGO, VIRGO and GEO600. Just recently, LIGO has 
reached its designed sensitivity and is currently undergo- 
ing enhancements to increase the sensitivity by an order 
of magnitude as a step towards advanced LIGO. In an- 
ticipation of these enhancements, it is essential to have 
models of gravitational waveforms for all sources of grav- 
itational radiation, in particular for binary black hole 
systems, since they are expected to be one of the most 
promising sources [H-Q . 

Constructing these waveforms is a non-trivial task. 
Generating a complete waveform involves numerically 
solving the full Einstein equations in order to correctly 
describe the last few orbits and merger. This is computa- 
tionally intensive with simulations running for weeks to 
produce a single accurate waveform. Furthermore, the 
parameter space of merging black hole binaries is quite 
large. In addition to the intrinsic black hole parame- 
ters (masses, spin magnitudes and orientations), there 
are the orbital parameters (eccentricity and semimajor 
axis). Because of the computational cost of producing 
numerical waveforms, the only way to have a hope of 
covering the parameter space efficiently is to use wave- 
forms that combine NR solutions with results from the 
PN approximation. To achieve this goal, it is firstly im- 
portant to cross-check the two methods to ensure that 
they give compatible results where PN is valid, namely 
for large enough binary separations. Secondly, it is nec- 
essary to investigate how close to the merger one can use 
the PN results (a recent study @ addresses this ques- 
tion in the extreme-mass-ratio case using the theory of 
optimal asymptotic expansions). 

Due to recent advances 0-Q hr the field of numerical 
relativity, long-term accurate and stable evolutions of bi- 
nary black hole systems spanning several orbits as well 
as the merger are now possible. The initial separations 
of the black holes in these simulations arc now sufficient 



to make comparisons with waveforms generated in the 
PN approximation. Such comparisons have been made 
for equal-mass non-spinning fll)| - fl7l ]. unequal- mass 
and spinning fl9| binaries, all in quasi- circular orbits. As 
a result, hybrid waveforms for quasi-circular orbits have 
been constructed which combine the PN waveform, ac- 
curate when the black holes are far apart, with the late 
inspiral and merger waveform that can only be obtained 
using full NR $M2^. 

In this paper, we take the next step and extend for the 
first time the NR and PN comparison to the case of eccen- 
tric binary black hole systems. Only recently the first NR 
studies of bound eccentric binary black hole systems have 
been performed 0, [25[ , where the dependence of the fi- 
nal black hole mass and spin on the initial eccentricity 
for non-spinning, equal-mass systems was studied. It has 
long been known that far-separated eccentric binary sys- 
tems emitting gravitational radiation will circularize [261 ] . 
However, it was not known what would happen if a bi- 
nary black hole system still had significant eccentricity 
in the late stages of inspiral. Rather than forming a final 
black hole with a higher or lower spin, the NR results 
showed that even up to an initial eccentricity of e ~ 0.4, 
the final black hole mass and spin were the same as in the 
circular case, indicating that the rate of loss of eccentric- 
ity was sufficient for the binary to circularize prior to or 
during merger. Due to the tendency of eccentric binary 
systems to circularize, most of the expected astrophysical 
binary black hole sources for Earth-based gravitational 
wave detectors will have lost all their eccentricity by the 
time their waves enter the frequency band of the detec- 
tor. However, several astrophysical scenarios have been 
proposed in which binary black hole systems in eccentric 
orbits might be detectable, for which it will be necessary 
to understand the dynamics and waveforms of eccentric 
binary black hole systems. 

One such scenario may occur in the dense cores of glob- 
ular clusters, where interactions between pairs of binary 
black hole systems eject one of the black holes, resulting 
in a stable hierarchical triple. This is a three-body system 
(three-body black hole systems have also been studied re- 
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cently in NR [27|) consisting of two closely bound black 
holes and a third orbiting the center of mass of the first 
two. When the two orbital planes are strongly tilted with 
respect to each other, tidal forces from the third body can 
cause an orbital resonance, increasing the eccentricity of 
the inner binary. This is known as the Kozai mecha- 
nism [28[ . It has been suggested (2|| that this could lead 
to eccentricities greater than about 0.1 at the time the bi- 
nary enters the frequency band of advanced ground-based 
detectors, followed by merger driven by gravitational ra- 
diation reaction. Stellar mass black hole binaries in glob- 
ular clusters are expected to have a thermal distribution 
of eccentricities [30[ , and intermediate mass black holes 
in globular clusters are expected to have eccentricities 
between 0.1 and 0.2 while they are in the frequency band 
of the LISA detector |3l( . Supcrmassivc black holes are 



also sources for LISA, and it is currently unknown what 
eccentricity they might have [32| • They could potentially 
merge within the Hubble time from hig hly eccentric or- 
bits if the Kozai effect was occurring [33[ . It has also been 
shown that massive black hole binaries in disks of gas can 
merge without losing eccentricity if the disc is rotating 
in the opposite sense to the binary orbit [34| • Being able 
to measure the non-zero eccentricity from the waveforms 
will tell us about the physics of the system, and may also 
have implications for detection if quasi-circular templates 
are used. 

The energy and angular momentum fluxes from the 
gravitational waves emitted by a comparable mass ec- 
centric binary were originally determined by Peters and 
Mathews [Iff, [35[ in the Newtonian limit. By balancing 
the time-averaged far-zone fluxes of energy and angular 
momentum with the loss of binding energy and angular 
momentum in the orbit, the rate of decay of the orbital 
scmimajor-axis and eccentricity could be determined in 
the adiabatic approximation. The result showed that the 
eccentricity of a binary reduces by approximately a factor 
of three when the semimajor axis is halved. 

The next order corrections to this result were obtained 
to 1 PN and 1.5 PN order, enabling the study of the evo- 
lution of the orbital elements using the quasi- Keplerian 
parametrization of the orbit [3a - l40j . With the use of a 
generalized Keplerian representation (41H43j . this work 
was extended to 2 PN ptl H^. 

An improved method of variation of constants has been 
developed [46|, |47} in order to construct models which 
for the first time go beyond the adiabatic approxima- 
tion. Very small oscillations in the orbital elements were 
found on the timescale of the orbital period. The con- 
servative 3 PN dynamics of an eccentric system in the 
quasi-Keplerian representation have been derived |48| . 
Recently, the complete 3 PN energy and angular mo- 
mentum fluxes have been determined p^ - IB"^ . 

The availability of the energy flux to 3.5 PN order [53[ 
in the quasi-circular case has led to successful matches 
with NR waveforms, with agreement in the waveform 
phase within 0.05 radians between 30 and 15 cycles before 
merger, and within several radians up to the merger 15 1 



for some PN models, though the level of agreement near 
merger is model dependent. The TaylorT4 model, specif- 
ically, agrees within 0.05 radians up to Muj gw = 0.1. Re- 
cently, it has been shown that for TaylorT4, the energy 
flux is identifiably different from the NR result even 25 
cycles before merger [l7| . 

A circular binary black hole inspiral gives rise to wave- 
form dynamics which are in some sense simple: the am- 
plitude and frequency increase monotonically, which may 
explain why the adiabatic approximation works so well. 
Eccentric orbits on the other hand give rise to waveforms 
with oscillations in the amplitude and frequency, and 
comparison with NR in this case will provide a signifi- 
cantly more stringent test of the PN approximation. 

In this paper, we present the first analysis of the agree- 
ment between PN and NR eccentric waveforms. We re- 
strict to the equal-mass, non-spinning case. We use the 
3 PN conservative quasi-Keplerian orbit equations |48j |. 
combined with the 2 PN evolution of the orbital ele- 
ments [46| to construct adiabatic PN waveforms, deter- 
mined by four independent initial parameters. We then 
present a full NR evolution which starts 21 gravitational 
wave cycles before merger with an estimated initial ec- 
centricity of e sa 0.1. We assume that the NR simulation 
gives the final stage of a full waveform, such as one that 
would be observed in nature. We then choose a fitting 
interval in time and use least squares fitting to find the 
parameters of the PN waveform which best matches the 
numerical data in that interval. We find agreement be- 
tween the NR and PN gravitational wave phase within 0.1 
radians for 10 wave cycles at the start of the simulation. 
The NR and PN phase difference grows to 0.6 radians 5 
cycles before merger, corresponding to Mlo. 
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As has been previously shown in the circular case, we 
find that different PN approximants lead to different lev- 
els of agreement with NR [l^]. We show here that an 
eccentric PN model expanded in terms of the mean mo- 
tion n = 2tt/P, where P is the orbital period, leads 
to significant disagreements with NR, whereas using the 
frequency-related variable x = ((2tt + A0)/P) 2 ^ 3 , where 
Acj) is the precession angle per period, gives much better 
agreement. 

In Section III A[ we describe the eccentric PN model 
we will use in our comparisons. The PN expressions are 
given in outline form to make clear precisely how we 
are constructing the solutions; the full expressions are 
given in the appendix. In Section III Bt we describe the 
methods used in our NR simulations that have not previ- 
ously been described; specifically, we discuss the method 
of constructing initial data parameters with eccentricity 
e pa 0.1. We present the results of the numerical simu- 
lations in Section 1111 A[ along with an analysis of the er- 
rors. In Section Hi CI we describe the method we use for 
matching NR and PN waveforms. Section [til Bl contains 
the main result of this paper, which is the comparison of 
the PN and NR solutions. Finally, we discuss the conse- 
quences of the results and our plans for future work in 
Section IT?] 
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II. METHODS 
A. Eccentric post-Newtonian model 

We first review the solution of eccentric Newtonian or- 
bits, in order to fix notation and to illustrate our gen- 
eral method for solving the PN system. For a detailed 
treatment, see for example Ref. [54] 0- The system under 
consideration consists of two point particles of masses mi 
and ru2- The total mass is M = m\ + m2- We will use M 
as the mass scale for all numerical quantities in our NR 
simulations, and work in units in which G = c = 1. The 
reduced mass is fx = m,\m-2.IM and the symmetric mass 
ratio is r\ = fx/M. We will give expressions for arbitrary 
mass ratios r], although in this work we will only be con- 
sidering equal-mass systems, for which mi = m% = M/2, 
li = M/4 and rj = 1/4. For Newtonian orbits, the energy 
E and angular momentum J are constants of the mo- 
tion and can be expressed in terms of the mean motion 
n and the eccentricity e. The conservation of J means 
that the orbit is restricted to a plane. The mean motion 
is related to the orbital (pericenter to pericenter) period 
P and the semimajor axis a by n = 2tt / P = aT^^M 1 / 2 . 
In the Newtonian case, the pericenter occurs at the same 
value of the relative angular coordinate 4> on each orbit; 
i.e. there is no precession. There is no closed form solu- 
tion for the relative orbital radius r or angular frequency 
(f> in terms of time, but they can be expressed in terms 
of the eccentric anomaly u, 



r = a [1 — e cos u] 
ny 



[1 — e cos u] 

The eccentric anomaly u satisfies Kepler's equation, 
I = u — e sin u , 



(1) 
(2) 

(3) 



where the mean anomaly I is given by I = n. Since n 
is a constant, we can integrate to obtain / = n(t — to) 
and Eq. ^ is a transcendental algebraic equation for 
u, which can be solved numerically, for example by New- 
ton's method, at each t. Thus we can obtain r and 4> (and 
hence f and <f>) at any time t. Each orbit is parametrized 
by the constants n, e, 4>o = 4>{fo) and Iq = l(to). 

The Newtonian system is conservative in the sense that 
it admits a conserved energy and angular momentum, 
which can be expressed in terms of the constants n and 
e. One can also derive conservative equations in the PN 
case; the Newtonian equations for r, <f> and I are modified 
by the addition of higher order (in n) terms. In the PN 
case, the quasi-Keplerian parametrization leads to three 
eccentricities, et, e r and e^, representing deviations from 



1 Note that in Ref. lo-ll . the eccentric anomaly is called tp rather 
than u, and the period T rather than P. The notation used in 
this work reflects that commonly used in the PN literature. 



circular motion in t, r and 4>, but these are related to each 
other by PN equations and it is sufficient to consider just 
e t - To Newtonian order, all three are equal. 

In the conservative PN equations, n and et are still 
constants, but the orbits precess. Note that the period 
P of the orbit is defined to be the time from pericenter 
to pericenter, and due to the effects of precession, this 
is not the time to go from angular coordinate <fi to <j) + 
2ir. The angle of precession of the pericenter during one 
(pericenter to pericenter) orbit of period P is denoted 
A<j>. Following Refs. 0,|H3], we define 



2?r + A(f> 
P 



(4) 



to be the angle swept out by the orbit from pericenter to 
pericenter in one period P. Note that in the conserva- 
tive PN system, this is a constant. In the circular case, 
where <j> is a constant, we have u> = 4>. We will investigate 
two different PN models which differ only in the choice 
of variable used (and hence in higher order uncontrolled 
remainder terms). In Ref. 543] , the eccentric system is 
described in terms of the mean motion n and the eccen- 
tricity et- We present the equations here in terms of the 
variable x = (Mw) 2 ^ 3 and et. We call the two resulting 
PN models the x-model and the n-model. See Sec. IIII CI 
for more discussion of these two models. 

We now give the 3 PN conservative orbital dynamics; 
we work throughout in modified harmonic coordinates, in 
which these expressions have been derived [48[ . The 3 PN 
conservative dynamics were first determined in Ref. [48j | , 
and were written out explicitly in terms of n and e t in 
Ref. [47|. Here, for brevity, we will omit lengthy high 
order PN expansions; the full expressions are available in 
the appendix. The abbreviated forms of the conservative 
dynamics, in terms of x and et, are 

r/M = [1 — et cos it] x~ x + ?tpn + f2PN# 

+ r 3PN x 2 + 0(x 3 ) (5) 



M(j> 



T x 3/2 + 0i PN a; 5/2 + 02PN£ 7/2 



[1 — et cos u] 1 

+ <fepNX 9/2 + 0(z 11/2 ) (6) 
I = u — et sin u + hp^x 2 

+ l 3PN x 3 + 0(x' 1 ) (7) 
Ml = Mn = x 3/2 + ni PN x 5/2 + n 2 pNX 7/2 

+ n 3PN x^ 2 + 0(x 11 ^ 2 ), (8) 

where the quantities tipn , </>ipn , • • ■ are functions of et 
and u, but tiipn, ■ ■ • arc functions only of et- Since the 
right hand side of Eq. ([5J is given in terms of the con- 
stants x and et, it can be trivially integrated to give l(t) 
in terms of an integration constant Iq at some to . So 
given constants x, et and Iq, we can solve Eq. ([7J nu- 
merically for u at each t by root-finding, then insert u 
into Eqs.( [5HS]) to obtain the coordinate motion of the 
conservative 3 PN system. 
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The conservative system is expected to be a good ap- 
proximation on timescales over which the energy and an- 
gular momentum lost to gravitational radiation is negli- 
gible. To go beyond this approximation, we will model 
these losses adiabatically; i.e. they will be averaged over 
the orbital period. The losses are derived by comput- 
ing the gravitational wave energy and angular momen- 
tum flux at infinity and equating the energy and angular 
momentum radiated to that lost from the system. The 
equations for E and J can be used to derive equations 
for x and e t . The equations to 2 PN order are given in 
Ref. |47| in terms of n and et. In terms of x and et, we 
have 



Mx 



Me 



2t? 



(96 + 292e 2 + 37e?) x 5 + x 1PN x 



' 15(l-e?)V2 
+ ii.5PN2: 13/2 + x 2PN x 7 + 0{x lb ' 2 ) (9) 

15(1^)5/2 ( 304 + 12l£ ?) + ^PNX 5 

+ ei. 5 PNX 11/2 + e 2PN z 6 + 0{x 13 ' 2 ) . (10) 



Since the evolution is treated adiabatically, the functions 
iiPN, cipn, • ■ • depend only on et, and not on u. Hence, 
the adiabatic evolution equations for x and e t form a 
closed system, and can be solved independently of the 
Kepler equation. Given initial conditions x(0) and et(0), 
we can solve the system of ODEs numerically to obtain 
x(t) and et(t). 

In the presence of time- varying x and et , Eq. ([8]) must 
be integrated to obtain l(t). The rest of the computation 
proceeds as in the conservative case; u is determined nu- 
merically by root- finding in Eq. ([7]), and then u, x and 
et are inserted at each time into Eqs.( [5] . 

We use analytical expressions for the functions 
'"lPN; 01PN) ■ • ■ to obtain numerical solutions for r and 
<fi, but due to the complexity of the expressions for r and 
0, we choose to obtain r and by numerically differenti- 
ating and integrating the r and solutions respectively 
This makes a difference to terms at higher PN orders that 
we are currently neglecting. 

We have checked our expressions for r and 0, as well 
as the 3 PN Kepler equation, by deriving them from the 
orbital elements in terms of E and h [48| , and comparing 
with the explicit expressions in terms of n and e t [47] ]. 
This completes the description of the coordinate motion. 

Since the NR and PN solutions are in different co- 
ordinate systems, we must compare them using some 
coordinate-independent quantity. We will use the grav- 
itational wave frequency; specifically the I = 2, m = 2 
mode of the Newman-Penrose ^4 quantity, as it is readily 
available from the NR simulations. 

The complex PN waveform strain is given (to leading 



Newtonian order) by 
h = h + — ih x 



(11) 



Mr} 



cos 2 + 1) 



-r 2 + r 2 2 + —\ 



+ 2rr0sin20' 
2Mr) 



• 2 2 12 

-r — r 



M 



sm 



■ cos ( 



R 

2r cos 2(f>' f 



M 



sin 2(j>' 



(12) 



(13) 



where 0' = — ip, and and ip are the spherical po- 
lar angles of the observer. Eqs. ITTUTBl are taken from 
Ref. [45( but with the sign convention for the cos 20' and 
sin 2^>' terms of Refs. 0, [47j. Using only the leading 
(quadrupolar) contribution to h is called the restricted 
waveform approximation. The strain h can be decom- 
posed into spin-weight s = —2 spherical harmonics, and 
the £ = 2, m — 2 mode is given by 



k = / -2Y2 (fl, <P)H0, <P)dM 




where 

coordinates, <p, tp, r, r, 
2, m = 2 spin-weight s 



(14) 
(15) 



/ 5/Vcos 4 (0/2). We insert the 
into Eq. (fi"5|) to obtain the I = 
= — 2 !53 mode of the waveform 



strain. Finally, using ^l 2 = h 22 we differentiate h 22 twice 



with respect to time to obtain the (complex) I = 2, m = 2 
mode of ^4. This is split into amplitude and phase, and 
undetermined additive multiples of 2-k in the phase are 
determined by continuity. 

We have described one procedure for constructing PN 
eccentric waveforms. Note that this is not unique; dif- 
ferent procedures will differ by the 'uncontrolled remain- 
der terms' of higher PN order than we have considered. 
Specifically, we have chosen to solve the 2 PN truncated 
adiabatic evolution equations for x and et numerically, 
rather than constructing an analytic expansion for the 
solution and then truncating it to 2 PN. This makes a 
difference to the solution in the circular case [ll| , and has 
been shown [l5[ to give better phase agreement with NR. 
In Ref. [IH , the circular waveform constructed using this 
approach is named TaylorT4, and the waveform phase 
agrees significantly better with NR than the TaylorTl, 
TaylorT2 and TaylorT3 approximants. For simplicity, 
we have also limited the computation of the waveform as 
a function of the coordinates to Newtonian (quadrupo- 
lar) accuracy, and restrict our comparisons to the phase 
rather than the amplitude. Higher order corrections to 
the waveforms are available [3 [56j], though not in a form 
which is convenient to use in this work. We have also cho- 
sen to construct some derivative quantities by numerical 
differentiation; where this is the case, we have verified 
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that the effects of discretization on the resulting wave- 
form phase are much smaller than any numerical errors 
we have in our NR simulations. 



B. Numerical relativity methods 

Our NR simulations are based on the moving punc- 
tures approach without excision 0, Q . Initial data rep- 
resenting the binary black hole system is constructed 
with a conformally flat metric and Bowen-York extrin- 
sic curvature, and the constraints are solved using the 
TwoPunctures [57j spectral code. The evolution in time 
is performed using our BSSN [5814601 finite differenc- 
ing code generated using the Kranc [61j code generation 
package. The Cactus [62j infrastructure is used for par- 
allclization, I/O and parameter hand ling , and for adap- 
tive mesh refinement we use Carpet [63(. The code has 
been previously described in more detail [64j], however 
we have since modified it to use sixth order spatial finite 
differencing as described in Ref. [65| in order to improve 
accuracy. We here use 9 levels of box-in-box mesh refine- 
ment, where the outermost (base) grid covers the domain 
x % G [—384, 384]. On the outer boundary, a simple spher- 
ical outgoing wave boundary condition is applied to each 
variable as is conventional for finite differencing BSSN 
codes in NR (see Ref. (6(| for more details). Formally, 
this boundary condition respects neither the constraints 
nor the characteristic structure of the equations. For 
very short simulations, it is possible to place the outer 
boundary far enough out that it is causally disconnected 
from the coordinate spheres on which the waveforms are 
computed, but for the long simulation we present here it 
is not computationally feasible to do this in our code. A 
discussion of the possible errors introduced can be found 
in Sec. ITiTAl 

The free parameters in the Bowen-York extrinsic cur- 
vature are the coordinate locations and linear momenta 
of the two black holes. We obtain these parameters us- 
ing the conservative 3 PN expressions for eccentric or- 
bits [43. These expressions require specification of the 
two constants, et and n (the eccentricity and mean mo- 
tion). We choose n = 0.0156/M and e t = 0.1, and 
compute the coordinate separation, r, from the 3 PN 
expression in terms of n and et, and use it in the Bowen- 
York extrinsic curvature. The tangential linear momen- 
tum of each black hole at apocenter, p y , is obtained from 
J = p y r, where J is computed as a PN expansion in n 
and e. We solve iteratively for the base mass parameters 
to ensure that the irreducible masses of the black holes 
are mi = m,2 = 0.5M, where M is a mass scale. As such, 
the mass scale M is the sum of the irreducible masses 
of the black holes. This procedure results in initial co- 
ordinate locations x l ± = (±7.1570737463, 0, 0)M, initial 
linear momenta P l ± = (0, ±0.07191137095, 0)M, and ini- 
tial bare masses m± arc = 0.4903157830M. The resulting 
spacetime has ADM mass A/adm = 0.991413M. 

This choice of initial data parameters has the follow- 



ing limitations. Firstly, only the conservative PN ex- 
pressions have been used, which means that there is no 
consideration of the inspiral velocity. Secondly, there will 
be an error in the parameters due to the truncation of 
the PN series. Thirdly, the use of PN parameters (in 
this case in harmonic coordinates) directly substituted 
into the Bowen-York extrinsic curvature, assumes that 
the differences in the coordinate systems are small. We 
will sec later that these initial data parameters agree rea- 
sonably well with the subsequent evolution. 



C. Fitting the post- Newtonian model to numerical 
relativity data 

We now discuss our method for determining a PN 
model which corresponds to our numerical simulation re- 
sults. The PN approximation is very accurate when the 
binary system is far separated, becomes less accurate in 
the later stages of inspiral, and is no longer valid during 
some period leading up to merger. Using NR, we can sim- 
ulate the late inspiral and merger. Ultimately, we would 
like to construct a waveform which most closely resembles 
one that would be observed in nature from early inspi- 
ral all the way through to merger. We will assume that 
the NR result gives the final part of this hypothetical full 
waveform, and use a PN waveform to approximate the 
full waveform before the start of the NR one. In this 
paper, we will not construct a hybrid waveform from the 
PN and NR results. 

In this work, we will look for agreement in the grav- 
itational wave frequency of the £ = 2, to = 2 mode of 

*4, 



4gw = dt arg * 



(16) 



as is common in the circular case. We will use the suffix 
'gw' to indicate that the quantity we are considering is 
related to the gravitational wave, and not the coordinate 
motion. We choose a time interval [t\, t%\ in the numerical 
simulation and use least squares fitting to determine the 
parameters of the PN model that best fits the numerical 
data in that interval. We will find in Section IIII Al that 
the black hole masses in the numerical simulations are 
essentially constant at mi = n%2 = 0.5M for the inspiral 
part of the simulation, so we do not fit for the masses 
when matching to PN. Thus, the eccentric PN model is 
determined uniquely by a choice of the functions X, et, 
I and 4> at a given time to, where X = x or n depending 
on the PN model being constructed. We define initial 
conditions 



Ho = [A o ,e o ^o,0o] 



(17) 



and the residual 



Q(yo) = ^2 [ wpN (^ y°) - w NR.(i)] 2 , (is) 



tei 
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using points t from the numerical simulation in the in- 
terval [tijta]- Q is then minimized numerically over y , 
where for each yo, the PN equations must be solved 
to construct the waveform. The minimization requires 
an initial estimate of yo- We find that using a local 
minimization method (for example, the principal axis 
method) can lead to inconsistent results. Specifically, 
the final fitted parameters show a dependence on the ini- 
tial estimate due to the existence of local minima in the 
residual. Instead, we use a global minimization method, 
requiring an order of magnitude more iterations (typi- 
cally around 5000), and hence increased computational 
resources. We find that minimization by the method of 
differential evolution, as implemented in Mathematica's 
NMinimize function, works well. A typical minimiza- 
tion for a given fitting window takes about 20 minutes 
on a laptop. Note that since the wave frequency w gw is 
independent of 4>o, in practice we determine <po by a sep- 
arate least squares fit between the PN and NR waveform 
phases. 

So, given a fitting interval, we can determine a PN 
model, identified by the parameters yo. If the model and 
data matched exactly, the fitted parameters yo would be 
independent of the fitting interval. However, the errors 
in the PN approximation cause the fit to be imperfect. 
If these errors are large, the dependence on the fitting 
interval will be significant. 

Once the parameters yo have been obtained, we can use 
these parameters to construct a final PN model, which 
will be the model that best approximates the full solution 
in the fitting interval. 



III. RESULTS 

A. Numerical relativity simulation results 

In this section, we describe the results of our NR simu- 
lations, and analyze the numerical errors. The PN model 
gives the limiting form of the waveform at large distances 
from the source, whereas in the numerical code we com- 
pute the waveform on coordinate spheres of finite radii 
r ext /M = {30, 40, ... , 150}. We therefore extrapolate the 
waveform to infinite radius using the method described 
in Ref. To extrapolate the waveform, we first shift 
the waveform measured at each extraction radius in time 
by the estimated light propagation time to the extrac- 
tion sphere, given by the Schwarzschild tortoise coordi- 
nate [671 ] . 

r* = r arcal + 2A/adm log ( J.T* 1 - , (19) 

\ ZM ADM J 

where we approximate r arca i w r + Madm (see Ref. [l5[ 
for further details; even if this relation does not hold 
exactly, the deviation will be included in our extrapola- 
tion error estimates). The amplitude and phase are then 
separately extrapolated by a least squares fit to an nth 



degree polynomial in 1/r, /™(r) = /oo + X)"=i a »/ ?,i at 
each time t — r*. We estimate the error in the nth order 
extrapolation as e" = f n+1 — f n . We find that using 
extraction radii r cxt = {70, 80, . . . , 150} in combination 
with first order extrapolation gives the best results. Us- 
ing higher order extrapolation does decrease the error, 
but the extrapolant contains more noise. 

We ran three simulations at different resolutions 
in order to assess the finite differencing error. The 
finest refinement boxes were coordinate cubes of side 
1.24M and consisted of 48 3 ,64 3 ,80 3 points in the three 
runs. This leads to finest grid spacings of hf = 
M/38.7,M/51.6,A//64.5. To investigate the finite dif- 
ferencing error, we consider the convergence properties 
of the gravitational wave phase, 

^(t -r*) = arg [*?(«- r*)] , (20) 

extrapolated to infinite radius. In Fig. [TJ we plot the 
convergence order of 4>. For t — r* < 500, we see no 
clear convergence order, but the differences between the 
phases at the three resolutions are less than 0.01 radi- 
ans. For 500 < t — r* < 2000 we see a convergence order 
which drops from 6 to 5, after which the order drops 
to about 1 for a small period around the merger. The 
fact that the convergence order is not clearly 6 may be 
explained by the fact that we have second, fourth and 
sixth order components in the simulation. Since we do 
not have clean sixth order convergence, we cannot reli- 
ably use Richardson extrapolation to obtain a more ac- 
curate result. However, we can use extrapolation of the 
highest two resolutions using the observed approximate 
convergence order of 5 to provide an estimate of the er- 
ror in the solution. Note that for the time region we will 
use for matching with PN (t < 1000), the convergence 
order of 5 is a good approximation. Figure [2] shows the 
finite differencing error estimate compared with the esti- 
mate of the error in the extrapolation to infinite radius. 
The dotted line represents the time of the peak in |\&4 2 |, 
which is a good indicator of the merger time. Note the 
sudden increase in the extrapolation error shortly after 
the merger. Also note that any significant effects arising 
from numerical reflections of the waves from refinement 
boundaries are expected to be covered by the finite dif- 
ferencing error bars, as these effects should diminish with 
increased resolution. 

When comparing with PN later, we will add the errors 
from finite differencing and from extrapolation in quadra- 
ture to provide an estimate of the overall error in the 
numerical waveform. Note that the approximately expo- 
nential growth of the finite differencing error in Fig . [5] 
has been previously observed in the circular case [65| . 

When comparing NR and PN models, we wish to 
use the gravitational wave frequency w™. However, as 
seen before for both finite differencing |l3| and pseudo- 
spectral [l5| codes in the circular case, ui gw has noticeable 
high frequency error at early times when the amplitude of 
the radiation is low. In our case, this comes from numer- 
ical reflections of the initial spurious radiation, present in 
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FIG. 1. Convergence order of the NR gravitational wave 
phase gw . Deviations from the expected value of 6 may be 
caused by lower order components in the code. 



FIG. 3. Filtering of NR gravitational wave frequency in the 
Fourier domain. The solution in the region containing noise 
is truncated to the lowest 30 Fourier modes. 
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FIG. 2. Errors in the NR gravitational wave phase </> gw from 
the effects of finite resolution and extrapolation to infinite 
radius 



the initial data, from mesh refinement boundaries. Since 
this is precisely the regime in which we would like to 
match with PN, this high frequency error must be re- 
moved. We find that this can be achieved very effec- 
tively by filtering the noisy region of ui gw in the Fourier 
domain. We first tried using a moving averages filter, 
but this tended to systematically reduce the amplitude 
of the oscillations in w gw , which we found unacceptable. 
We also chose not to fit a polynomial to w gw as has been 
done in the circular case due to the naturally os- 
cillatory nature of the eccentric signal. To perform the 
filtering, we proceeded as follows. We first chose an in- 
terval of time, [^1,^2], m which to perform the filtering. 
We chose [h,h] = [80 M, 1680 A/], excluding the late in- 
spiral and merger as well as the initial spurious radiation 
from the filtering region. We then performed a discrete 



Fourier transform of the data, removed all but the low- 
est 30 modes, and then inverse transformed. We found 
that 30 modes were sufficient to represent the signal; this 
was judged by subtracting the filtered from the unfiltered 
signal, and observing essentially only noise. Taking only 
the first 30 modes corresponds to a frequency cutoff of 
u mal = 30 x 2-K/{t 2 — t\) w 0.1M -1 , or modes with a 
period of T m - m ps 50M. Note that this is not comparable 
to filtering the evolved variables or even vE^ 2 i it is the 
frequency of ^\ 2 that is being filtered. Since the original 
signal is not periodic, Gibbs phenomena were observed as 
oscillations near the endpoints of the filtered region. We 
therefore removed a segment of length 80M from the be- 
ginning and end of the filtered region before re-inserting 
the filtered region into the full signal. Fig. [3] shows the 
result of the filtering. 

We have monitored the irreducible masses M- m of the 
apparent horizons in the lowest resolution simulation. 
The computed mass of each black hole drops from its 
initial value of 0.5 by only 2 x 10~ 4 M by the time of the 
merger, and we ascribe this effect to finite differencing 
error. We have not computed horizons at higher reso- 
lutions due to computational expense. Thus, within our 
numerical errors, we do not detect any physical growth of 
the horizons during the inspiral, which potentially could 
have occurred due to absorption of gravitational wave 
energy in the initial part of the simulation, as has been 
studied in detail in previous work (68|. 

The spins of the black holes, as measured using an 
approximate technique derived from the isolated horizon 
formalism [69l.[70l|. increase during the simulation to only 
S z = 10~ 4 /M 2 before the merger. This is independent 
of finite differencing resolution, but we expect this tiny 
spin to be of little consequence to the PN comparison, 
which does not contain the effects of spin. 

The outer boundary in the simulations is at x l = 
±384M, and as mentioned in Sec lII B[ the boundary con- 
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dition is a source of error in the simulation. To measure 
the effect of this error, we have repeated the low reso- 
lution simulation, which has only modest computational 
cost, with the outer boundary moved to x l = ±768 A/ 
by enlarging the coarsest grid. We find that the effect 
on the waveform phase is much smaller than the esti- 
mated errors in the high resolution simulation due to fi- 
nite differencing and extrapolation to infinite radius, and 
we conclude that the outer boundary is not a significant 
source of error in the simulation. In future, with more 
accurate simulations, this will need to be addressed fur- 
ther. 

The simulations at the three different resolutions con- 
sumed approximately 5000, 11000 and 16000 CPU hours 
respectively, each one running on 32 cores of the LoncStar 
supercomputer. 



B. Comparing numerical relativity simulations 
with post-Newtonian models 

We now discuss the results of applying the fitting pro- 
cedure described in Section III CI to the numerical simula- 
tion results. 

Figures 0] and [5] show the parameters [xq, eo, h, <f>o] (f° r 
the x model) and [no, eo, lo, <f>o] (for the n model) deter- 



mined by fits of the NR data to the PN model in fitting 
intervals / = [ii,^]- These parameters are the values of 
the functions x, n, e, I and <f> at t — r* = 0. In Fig. [U t\ 
has been kept fixed to a value at the start of the usable 
waveform and £2 has been varied. We see that the pa- 
rameters obtained from fits using the x model vary much 
less with the fitting window length than those using the 
n model. Specifically, we see that for both models the fit- 
ted parameters oscillate significantly for intervals of less 
than ~ 400M, but for the x model these variations die 
away as the interval is increased beyond this. From the 
initial data parameters, the orbital period is P — 403Af. 
It may be that over timescales smaller then the orbital 
period, there are unmodeled non-adiabatic oscillations in 
the NR result which are averaged out when larger fitting 
intervals are used. These oscillations may cause the fit 
to become worse for small intervals. For the n model 
we see strong oscillations of a period ~ 400A/ roughly 
corresponding to the period of the oscillations in uj gw it- 
self. In order to determine the effect on the parameters 
of the interval location, we choose an interval width of 
400Af and vary t\ in Fig. [5] Here again we see that the 
x model shows much more consistent behavior than the 
n model. 

In order to choose a unique set of PN parameters, we 
choose the earliest possible fitting interval, and take the 
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Parameter 


a;-model fit 


n-model fit 


Initial data value 


x 


0.0747729 


- 


0.0740853 


n 




0.0158959 


0.0156 


eo 


0.103291 


0.10299 


0.1 


lo 


3.06358 


2.9529 


7T = 3.1416 


<f>0 


-1.47386 


-1.45652 






TABLE I. Eccentric PN (a;-model and n-model) parameters 
computed by fitting in an interval [210, 610] as well as the 
parameters estimated from the initial data. The parameters 
correspond to the values of the functions x, n, e, I and <j> at 
t — r* = 0. Note that the agreement is not expected to be 
exact . 



size of the interval to approximately correspond to the 
initial orbital period, ~ 400M, giving a fitting interval 
t/M £ [210,610]. The parameters for this fitting inter- 
val are given in Table HI It is interesting to compare 
these parameters with the approximate parameters used 
to construct the Bowen-York initial data; these are also 
given in the table, xq and no agree to within 1% and 2% 
respectively with the initial data values, eo agrees within 
0.3% between the two PN models, and to 3% with the 
initial data value. Iq agrees to within 0.1 radians between 
the two models and the initial data value. 0o agrees to 
within 0.02 radians between the two models, but is of the 
order of 7r/2 different from the initial data value. This 
large discrepancy is probably related to the adjustment 
of the coordinate system that happens at the start of the 
numerical simulation. Recall that the method for con- 
structing the initial data parameters was approximate, 
due to the different coordinate systems used, so perfect 
agreement is not expected. 

Now that we have estimated the PN model which 
matches the NR solution in the fitting interval, we can 
compare the PN waveform for the x and n models with 
the NR result. In Fig. [BJ we plot the PN and NR gravi- 
tational wave frequencies co sw and see that there is good 
agreement with the x model from the start of the simu- 
lation to t w 1800M. That there is such a high level of 
agreement with a model which contains so much struc- 
ture is a strong validation of both the PN model and the 
NR simulation. We also see on the same plot the much 
worse agreement obtained using the n-model. 

We now quantify the agreement with the x and n- 
models by considering the waveform phase differences. 
Fig.[7]shows the difference between the NR and PN grav- 
itational wave phases as a function of t. The error bars 
represent the uncertainty in the NR phase from extrap- 
olation to infinite radius and finite differencing trunca- 
tion error. We see that the phase difference between NR 
and PN is within 0.1 radians for approximately 1330A-/, 
or 11 GW cycles. At t = 1882M, corresponding to 
Mw gw = 0.1, the phase difference between NR and PN 
is w 0.7 radians. 

To put the phase difference of 0.7 radians at Muj gw = 
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the n-model is significantly worse. 
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NR simulation and the PN a;-model. The error bars represent 
the estimated errors in the NR simulation. 



0.1 into context, we note that the TaylorT4 circular PN 
model, which is very similar to our eccentric model with 
e = 0, has been shown to have a phase difference at 
Mw gw = 0.1 of ~ 0.3 radians for 2 PN radiation reaction 
(see Fig. 22 in Ref. [HI). We should be cautious about 
drawing the conclusion that the agreement in the circular 
case is better, however, as AIu} gw = 0.1 may not be di- 
rectly comparable in the two cases, particularly because 
o; gw oscillates in the eccentric case, but is monotonic in 
the circular case. The steepness of the phase difference 
in Fig. 22 in Ref. [ll| at that point makes the comparison 
very sensitive to the exact point chosen. 
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C. Choice of post-Newtonian variables 

Throughout this work we have presented the results of 
fitting two PN models with NR data. The two models 
differ only in the choice of variable used: the frequency- 
related variable x or the mean motion n. Our first at- 
tempts at matching the NR simulation with an eccentric 
PN model used n. We studied this case extensively, but 
found significant disagreement, as has been shown. Faced 
with this disagreement, we studied the (much simpler) 
circular case using a simulation |25| with low-eccentricity 
initial data [7l[ and a circular PN model formed by tak- 
ing our eccentric n-model and setting e = 0. This model 
is suboptimal as it only has 2 PN radiation reaction, and 
3.5 PN expressions are available for the circular case. The 
agreement between NR and PN is very poor even in the 
circular case using n; the gravitational wave phase dif- 
ference at Muj gw = 0.1 is ~ 20 radians. (Note that one 
should be careful about making direct detailed compar- 
isons between the circular and eccentric cases, due to the 
ambiguity in the choice of reference point AIuj gw = 0.1 
due to the eccentric oscillations in w gw .) However, ex- 
pressing the PN equations in terms of the coordinate an- 
gular velocity of the black holes, w, as is common in the 
literature, gives a significant improvement over using n; 
at Mw gw =0.1, the phase difference is 0.8 radians. This 
is in broad agreement with the difference of ~ 0.3 radians 
in Fig. 22 of Ref. [H for the TaylorT4 model at 2 PN, ac- 
counting for the uncertainty in the choice of comparison 
time. This motivated us to search for a frequency-related 
variable applicable in the eccentric case, and we chose to 
use x = (Mw) 2 ' 3 . for compatibility with Ref. [4|| (recall 
that in the eccentric case, uj = (2ir + A(j))/P ^ <ft), lead- 
ing to the 0.7 radian phase difference at Muj gw =0.1 we 
report here. 



of the same order have been shown to have different errors 
in the NR regime. In particular, the TaylorT4 circular 
model showed a remarkable agreement in the waveform 
phase, but there was a noticeable disagreement in the 
energy flux [l7|. It has also been shown that this re- 
markable a gre ement is lost when spinning systems are 
considered [13| . Our eccentric PN model based on x is 
very similar to TaylorT4 as e — > 0, so we would expect 
the same conclusions to apply. 

Now that it is possible to match NR and PN eccentric 
waveforms, we plan to start to construct hybrid tem- 
plates and begin to assess the implications for the inter- 
ferometric detection of gravitational wave signals from 
eccentric binaries close to and including merger. Since 
complete 3 PN radiation reaction terms for the angular 
momentum flux have now also been computed, we will be 
able to compare with a fully 3 PN model, and expect the 
agreement with NR to get better closer to the merger. 
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IV. CONCLUSIONS 



We have presented NR results for an inspiraling ec- 
centric black hole binary system with initial eccentricity 
e « 0.1 and compared them with two adiabatic eccen- 
tric PN models {x and n) with 2 PN radiation reaction. 
For the x model, the gravitational wave phase agrees 
to within ±0.1 radians between 21 and 11 cycles before 
merger. The difference grows to 0.7 radians at ~ 5 cycles 
before merger (Mw gw = 0.1), in broad agreement with 
the circular case at 2 PN order. One cycle before the 
merger, the solution to the PN ODEs diverges, indicat- 
ing a breakdown of the model. 

We found that it was necessary to express the PN 
model in terms of the frequency-related variable x rather 
than the mean motion n to get this level of agreement. 
We conjecture that, when expressed in terms of n, cer- 
tain higher order PN terms are non-negligible, whereas 
when expressed in terms of a;, they are small, leading to 
a smaller error in the PN solution. This can be likened to 
studies [HI, [l5| where different circular PN approximants 
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Appendix A: PN expressions 



We now present, for reference, the full PN expressions used in this work. The expressions for the 3 PN conservative 
dynamics (i.e. r, <fi, Z, n) can be derived in two ways from the existing literature. They are given directly in Ref. (47j in 

terms of n and et, so all that remains is to express them in terms of x and e t . Recall that x is defined as x = (Mw) 2 ' 3 
where cj = (2tt + A(j))/P and P = 2ir jn. In Ref. [48| . $ is used in place of A<p, where $ = 2n + A<j>. This reference 
gives expressions for n, e t and $ in terms of E and J; these can be used to obtain n in terms of x and et, 
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where, for brevity, we have written e = e t - This expression for n is then substituted into the conservative expressions 
in Ref. 47J to obtain the conservative expressions in terms of x and et, dropping any resulting terms which are 
higher than 3 PN. Alternatively, we can derive these expressions by taking the expressions for the orbital elements 
in Ref. (48|, along with the expressions for r and (j>, all in terms of E and J. By both methods, we obtain for the 
separation r, 
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The relative angular velocity <fi is found to be 
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113680?; - 221760)e 9 + ( - 11200?; 3 - 112000?; 2 + 12915tt 2 ?; + 692928?; - 194880)e 7 + (4480?? 3 + 8960?; 2 
43050tt 2 ?; + 1127280?;- 147840)e 5 ) cos 5 (u) + (( - 16240?? 3 + 12880?; 2 + 18480??)e 10 + (16240?; 3 - 91840?? 2 
172207T 2 ?; - 652192?; + 100800)e 8 + ( - 55440?; 3 + 34160?; 2 - 30135tt 2 ?7 - 2185040?? + 2493120)e 6 
(21840?; 3 + 86800?; 2 + 163590tt 2 ?; - 5713888?; + 228480)e 4 ) cos 4 (it) + ((560?; 3 - 137200?; 2 + 388640?; 
241920)e 9 + (30800?? 3 - 264880?; 2 - 68880tt 2 ?; + 624128?; + 766080)e 7 + (66640?; 3 + 612080?; 2 - 8610tt 2 ?? 
6666080?; - 6652800)e 5 + ( - 30800?; 3 - 294000?; 2 - 223860tt 2 ?; + 9386432??)e 3 ) cos 3 (it) + 67200?; 2 
((4480?? 3 - 20160?] 2 + 16800?;)e 10 + (3920?; 3 + 475440?? 2 - 17220tt 2 ?; + 831952?? - 725760)e 8 + ( - 75600?; 3 
96880?; 2 + 154980tt 2 ?; - 3249488?; - 685440)e 6 + (5040?; 3 - 659120?; 2 + 25830tt 2 ?; - 7356624?; - 



( - 5040?] 3 
424480?;) e 1 



190960?; 2 - 
- (28560?; 3 



1377607T 2 ?; - 7307920?; + 107520)e 2 ) cos 2 (?i) - 761600?;+ (( - 2240?; 3 - 168000?; 2 
242480?; 2 + 34440tt 2 ?; - 1340224?; + 725760) e 7 + ( - 33040?; 3 - 754880?; 2 



1722007T 2 ?; + 5458480?; - 221760)e 5 + 



738640?? 2 + 301357T 2 ?; + 1554048?; - 2936640)e 3 



( - 560?? 3 - 100240?; 2 - 43050tt 2 ?; + 3284816?; - 389760)e) cos(it) + \/l -e 2 ((( - 127680?; 2 + 544320?; 
739200)e 7 + ( - 53760?; 2 - 8610tt 2 ?? + 674240?; - 67200)e 5 ) cos 5 (u) + ((161280?; 2 - 477120?; + 537600)e s 
(477120?; 2 + 17220tt 2 ?; - 2894080?; + 2217600)e 6 + (268800?; 2 + 25830tt 2 ?; - 2721600?; 
1276800)e 4 ) cos 4 (w) + (( - 524160?; 2 + 1122240?; - 940800)e 7 + ( - 873600?; 2 - 68880tt 2 ?? + 7705600?; 



3897600)e t 



416640?; 2 - 17220tt 2 ?; + 3357760?; - 3225600)e 3 ) cos 3 (i?) + ((604800?; 2 - 504000?; 



403200)e 6 + (1034880?; 2 + 103320tt 2 ?; - 11195520?; + 5779200)e 4 + (174720?; 2 - 17220tt 2 ?? - 486080?; 
2688000)e 2 ) cos 2 (it) + (( - 282240?; 2 - 450240?; + 1478400)e 5 + ( - 719040?; 2 - 68880tt 2 ?; + 8128960?; 
5040000)e 3 + (94080?; 2 + 25830tt 2 ?; - 1585920?; - 470400)e) cos(u) - 67200?; 2 + 761600?; + e 4 (40320?? 2 
2 (208320?; 2 + 17220tt 2 ?; - 2289280?; + 1680000) - 8610??7r 2 - 201600) + 8610?;7r 2 



309120?; - 672000) 
201600 . 



(A14) 
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The 3 PN Kepler equation is 

I = 'opn + ^2pn£ 2 + kpNX 3 + 0(x 4 ) 
?opn = u — e sin u 

1 



?2PN = 
fepN = 



8V1 - e 2 (l - ecos(w)) 
1 



[— 12(2?? — 5)(w — w)(ecos(u) — 1) — e\J 1 — e 2 (?? — 15)?? sin(it)j 

35(96(11t7 2 - 2977 + 30)e 2 + 960?? 2 +T}(- 13184 + 123tt 2 ) 



(A15) 
(A16) 

(A17) 



6720(1 - e 2 ) 3 / 2 (l - ecos(w)) 3 
+ 8640) (u - u)(ecos(u) - l) :i + 3360( - 12(2r? - 5)(u- «) + 12e(2?? - 5) cos(m)(m - v) 



+ ey/l - e 2 (r/ - 15)77 sin(w)) (e cos(u) - l) 2 + ei/l - e 2 (l40(l3e 4 - lie 2 - 2)t? 3 - 140(73e 4 - 325e 2 + 444)?? ; 
+ (3220e 4 - 148960e 2 - 4305tt 2 + 143868) 77 + e 2 (l820(e 2 - l)?? 3 - 140(83e 2 + 109)?? 2 - (ll20e 2 + 4305tt 2 
+ 752)?7 + 67200) cos 2 (tj) - 2e(l960(e 2 - l)rf + 6720 (e 2 - 5)?? 2 + ( - 71820e 2 - 4305tt 2 + 69948) 77 

+ 67200) cos(tj) + 67200) sin(i 
where, as in Ref. [13], we use 

sin (it 



v — u = 2 tan 1 



1 — cos(t 



(A18) 



(A19) 



and 



l-Jl-e3 



e<t, is given by 

= e + e^ipN^ + e^PNX 2 + e^3PNi 3 + 0(x 4 ) 
e^iPN = —e(r) - 4) 
e 



e</,2PN — 



96 e 2 - 1 



e</)3PN — 



(4I77 2 - 65977 + 1152) e 2 + 4?? 2 + 6877 + \/l - e 2 (28877 - 720) - 1248 

(1344077 2 + 4838407? - 940800) e 4 + (255360r? 2 + 172207r 2 7? - 28806407? + 2688000) e 2 



(A20) 



(A21) 
(A22) 

(A23) 



S/2 



26880(1 - e 2 ) 

- 2688OO77 2 + 239680077 + \]\ - e 2 ((l050?7 3 - 13405077 2 + 78631077 - 860160) e 4 + ( - I89OO77 3 + 55398077 2 
+ 4305vr 2 77 - 1246368?7 + 2042880) e 2 + 276640t7 2 + 2674480?? - 172207?7r 2 - 1451520) - 17220t?7t 2 

- 1747200 . (A24) 

This completes the expressions used in the conservative dynamics. The radiation reaction is given to 2 PN order in 
Ref. [47l | in terms of n and e t . We again substitute for n in terms of x, and obtain 



Mi = ioPNZ 5 + iiPN^ 6 + ±i.5Pn:e 13/ ' 2 + ±2pn£ 7 + 0(x 15/2 ) 
2(37e 4 + 292e 2 + 96)t? 



XQPK = 



(A25) 
(A26) 



15(l-e 2 ) 7/2 

*ipn = 77^ V „ a/ , \-(S2S8tj - 11717)e 6 - 14(10122?? - 12217)e 4 - 120(1330?? - 731)e 2 - 16(924?? + 743)1 (A27) 
420(1 — e 2 ) M ' 2 

256 

xi.spn = -=-7?vrK E (e) (A28) 



i 2PN = 45360(2 _ e 2)ii/2 (!964256?? 2 - 3259980?? + 3523113) e 8 + (64828848?? 2 - 123108426?? + 83424402) e 
+ (16650606060?? 2 - 207204264?? + 783768) e 4 + (61282032?? 2 + 15464736?? - 92846560) e 2 + 1903104?? 2 
+ i/l - e 2 ( (2646000 - 1058400??)e 6 + (64532160 - 25812864??)e 2 - 580608?? + 1451520) + 4514976?? 

- 360224 , 



(A29) 



14 



for x, and 



,,. • 4 , ■ 5 , ■ 11/2 . ■ 6 i 13/2\ 

Me = eopNX- + eipNX- + ei.sPNa: + e2PN2: + 0(x ' ) 
e(l21e 2 + 304)?? 



6opn 



eiPN = 



15(1 - e 2 
erj 



5/2 



ei.5PN 



2520(1 - e 2 
128??7r " 



7/2 



(93184?7 - 125361)e 4 + 12(54271?? - 59834)e 2 + 8(28588?? + 8451) 



5e 



(e 2 - l)Ks(e) + \/l-e 2 K j(e) 



(A30) 
(A31) 

(A32) 

(A33) 



e2PN = 



07/ 



n/2 



(2758560?? 2 - 4344852?? + 3786543) e 6 + (42810096?? 2 - 78112266?? + 46579718) e 4 



30240(1 - e 2 ) 

+ (48711348?? 2 - 35583228?? - 36993396) e 2 + 4548096?? 2 + \/l - e 2 ((2847600 - 1139040??)e 4 + (35093520 
- 14037408??)e 2 - 5386752?? + 13466880) + 13509360?? - 15198032 



(A34) 



for e. These equations are written in terms of the functions ke and kj, given in Ref. [46[ in terms of infinite sums 
of Bcsscl functions. We reproduce them here for completeness. 

KE = £ JP 3 ((( - e ' - i + i + 3 )P 2 + \ ~ 3 + i)Mpe) 2 + ( - 3e - — + 7 -)pJ' v {pe)J p {pe) 



P =i 



+ (( e 2 + ^-2)p 2 + 7>-i)^) 2 ) 



*J = E \v V / T^(( - -J - 1 + %)vUve? + (2(e + 1 - V - I + ^) Jj(pe) J p (pe) + 2(1 - ^)K(pe 



(A35) 



(A36) 



P =i 



These are functions of e only, and are computed numerically using a sufficient number of terms in the summation 
that the result converges to within machine precision (10~ 15 ). For computational efficiency, the resulting function 
is converted into an interpolating polynomial, and the interpolation error is estimated to be ~ 10 -12 in the range 
< e < 0.4. 
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